ℹ Using "','" as decimal and "'.'" as grouping mark. Use `read_delim()` for more control.

1 Etape 1 : Première manipulation du RPG

1.1 Je choisis un point sur la carte : https://www.google.fr/maps

=> clic droit sur la carte + clic gauche sur les coordonnées (copiées dans le presse-papier)

J’indique un rayon en mètres (< 15000) pour sélectionner les parcelles autour du point

Code
# coord_gmaps<- rstudioapi::showPrompt(title = "Collez les coordonnées", message = "coordonnées Gmaps", default = "")
coord_gmaps<-"46.46946179131805, -1.4932577154775046"
lat<-as.numeric(str_split(coord_gmaps,fixed(","),simplify = TRUE)[,1])
lon<-as.numeric(str_split(coord_gmaps,fixed(","),simplify = TRUE)[,2])
# rayon<-as.numeric(rstudioapi::showPrompt(title = "Rayon", message = "Rayon (en m)", default = ""))
rayon<-5000

# création d'une table sf «point» avec les coordonnées saisies
# transformation des coordonnées en syst de proj 2154 (Lambert II - Français) 
point<-data.frame(lon,lat,rayon) %>% 
  st_as_sf(coords = c("lon","lat"),crs = "EPSG:4326") %>%
  mutate(coord_pt_gps=st_as_text(geometry)) %>% 
  st_transform("EPSG:2154") %>% 
  st_sf() %>% clean_names() %>% 
  rename(geom=geometry)

# st_crs(point)
# st_geometry_type(point)
# plot(point)

La table des parcelles agricoles 2021 se trouve sur un serveur PostgreSQL/PostGIS.

1.2 Je me connecte au serveur PostgreSQL avec le mot de passe disponible sur Onyxia/Mes services/PostgreSQL/Read me

Code
# le mot de passe est stocké dans un secret Vault

# postgresql_password <- rstudioapi::askForPassword(prompt = "Entrez le password PostgreSQL")
postgresql_password <- "1tfawt3nj7fgzo3w7cma"
  
# Connection à PostgreSQL
cnx <- dbConnect(PostgreSQL(),
                user = "projet-funathon",
                password = postgresql_password,
                host = "postgresql-438832",
                dbname = "defaultdb",
                port = 5432,
                options="-c search_path=rpg,public") # specify what schema to connect to

1.3 je lance une requête SQL pour sélectionner les parcelles situées autour de mon point dans le rayon choisi.

Code
# suppression de la table «point» si elle existe
dbSendQuery(cnx,"DROP TABLE IF EXISTS rpg.point CASCADE;")
<PostgreSQLResult>
Code
# écriture de la table point dans la base PostGis
write_sf(point, cnx, append = T)

# ajout d'une clé primaire
dbSendQuery(cnx,"ALTER TABLE rpg.point ADD CONSTRAINT point_pkey PRIMARY KEY(coord_pt_gps);")
<PostgreSQLResult>
Code
# ajout d'un index
dbSendQuery(cnx,"CREATE INDEX ON rpg.point USING gist(geom);")
<PostgreSQLResult>
Code
# dbExecute(cnx,"CREATE INDEX ON rpg.point USING gist(geom);")

# 3 - Exécution de la requête de découpage du RPG autour du point sur PostGis  -------

dbSendQuery(cnx,"DROP TABLE IF EXISTS rpg.parc_prox CASCADE;")
<PostgreSQLResult>
Code
dbSendQuery(cnx,"CREATE TABLE rpg.parc_prox AS
    SELECT row_number() OVER () AS row_id, p.coord_pt_gps, p.rayon, r.*  
    FROM rpg.point p, rpg.parcelles r 
    WHERE ST_DWithin(p.geom,r.geom,p.rayon);")
<PostgreSQLResult>
Code
# ajout d'une clé primaire
dbSendQuery(cnx,"ALTER TABLE rpg.parc_prox ADD CONSTRAINT parc_prox_pk PRIMARY KEY(id_parcel);")
<PostgreSQLResult>
Code
# ajout d'un index
dbSendQuery(cnx,"CREATE INDEX parc_prox_geom_idx ON rpg.parc_prox USING gist(geom);")
<PostgreSQLResult>
Code
# 4 - Téléchargement des parcelles proches sous R-------------------------------

# parc_prox<-read_sf(cnx,parc_prox)

parc_prox<-st_read(cnx, query="select * from rpg.parc_prox;")

1.4 J’affiche les parcelles autour de mon point avec une carte interactive leaflet

Code
#plot(st_geometry(parc_prox))
#plot(st_geometry(point),add = T,col = "red")

# Transformation de la projection car leaflet ne connait que le WGS 84
carte_parc_prox<- parc_prox %>% st_transform(4326)

# Marqueur du point
pt_mark<- point %>% st_transform(4326)

# ajout du libellé des cultures
carte_parc_prox_lib<-carte_parc_prox %>% 
  left_join(lib_cult %>% select(-code_groupe_culture),by=c("code_cultu"="code_culture")) 
#création d'un label ad hoc à afficher en surbrillance au passage de la souris sur la carte.
labels<-sprintf("<strong>id_parcel : </strong>%s<br/>
                <strong>Groupe culture : </strong>%s<br/>
                <strong>Culture : </strong>%s<br/>
                <strong>Surface (ha) : </strong>%s<br/>
                <strong>Département : </strong>%s<br/>
                <strong>Commune : </strong>%s<br/>",
                parc_prox$id_parcel,carte_parc_prox_lib$libelle_groupe_culture,
                carte_parc_prox_lib$libelle_culture,parc_prox$surf_parc,
                parc_prox$insee_dep,parc_prox$nom_com) %>% 
  lapply(htmltools::HTML)

# labels

# création d'une palette de couleurs associée au groupe de culture
factpal <- colorFactor("Paired", parc_prox$code_group)

carte_parc_prox_html <- leaflet(carte_parc_prox_lib) %>% 
   # addProviderTiles("Esri.WorldImagery") %>%
   # addTiles() %>% 
  addTiles("http://wxs.ign.fr/essentiels/wmts?REQUEST=GetTile&SERVICE=WMTS&VERSION=1.0.0&STYLE=normal&TILEMATRIXSET=PM&FORMAT=image/jpeg&LAYER=ORTHOIMAGERY.ORTHOPHOTOS&TILEMATRIX={z}&TILEROW={y}&TILECOL={x}") %>%
  addPolygons( #fillColor="white",
               fillColor=~factpal(code_group),
               weight = 2,
               opacity = 1,
               color = "#ffd966",
               dashArray = "3",
               fillOpacity = 0.5,
               highlight = highlightOptions(
                 weight = 5,
                 color = "#A40000",
                 dashArray = "",
                 fillOpacity = 0.0,
                 bringToFront = TRUE),
               label = labels,
               labelOptions = labelOptions(
                 style = list("font-weight" = "normal", padding = "3px 8px"),
                 textsize = "15px",
                 direction = "auto",
                 encoding="UTF-8")) %>% 
  addMarkers(data=pt_mark,~lon,~lat, popup = ~coord_pt_gps, label = ~coord_pt_gps)
Warning in RColorBrewer::brewer.pal(max(3, n), palette): n too large, allowed maximum for palette Paired is 12
Returning the palette you asked for with that many colors
Code
# leaflet(pt_mark) %>% addTiles() %>% addMarkers(~lon,~lat)

carte_parc_prox_html
Code
# saveWidget(widget = carte_parc_prox_html, file = "carte_parc_prox_html.html")

1.5 Quelle est la structure des parcelles agricoles autour de mon point ?

Code
t1 <- parc_prox %>% st_drop_geometry() %>% count(code_group) %>% 
  add_tally(n) %>% 
  mutate(n_pct=round(100*n/nn,1)) %>% 
  select(-nn) %>% rename(n_parcelles=n) %>%
  # adorn_totals() %>% 
cbind(
  # comptage des surfaces
parc_prox %>% st_drop_geometry() %>% count(code_group,wt=surf_parc) %>% 
      add_tally(n) %>% 
      mutate(surf_pct=round(100*n/nn,1)) %>%
      select(-nn) %>%  
      rename(surf_parc_ha=n) %>% select(surf_parc_ha, surf_pct) # %>% 
       # adorn_totals()
) %>% left_join(lib_group_cult,by=c("code_group"="code_groupe_culture")) %>% 
  select(code_group,libelle_groupe_culture,everything()) %>% 
  # arrange(as.numeric(code_group)) %>% 
  arrange(desc(surf_parc_ha)) %>% 
  adorn_totals() %>% 
  mutate(taille_moy_parc=round(surf_parc_ha/n_parcelles,1))
Storing counts in `nn`, as `n` already present in input
Storing counts in `nn`, as `n` already present in input
ℹ Use `name = "new_name"` to pick a new name.
Code
t1  %>% 
  setNames(c("code","groupe de cultures","nombre de parcelles","(%)","surface (ha)","surface (%)","taille moyenne (ha)")) %>% 
  kable(
    format="html",
    caption="<span style='font-size:medium'>Groupes de cultures <strong>locales</strong> par surfaces décroissantes</span>",
    format.args = list(decimal.mark = ",", big.mark = " "),
    booktabs = TRUE) %>%
  kable_styling(font_size = 15) %>% 
  gsub("font-size: initial !important;",
       "font-size: 20pt !important;",.)%>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>% 
  row_spec(nrow(t1), bold = T, color = "white", background = "grey")
Structure des cultures au niveau local
code groupe de cultures nombre de parcelles (%) surface (ha) surface (%) taille moyenne (ha)
2 Maïs grain et ensilage 306 18,8 1 669,06 30,4 5,5
18 Prairies permanentes 402 24,7 1 168,78 21,3 2,9
1 Blé tendre 178 10,9 737,43 13,4 4,1
19 Prairies temporaires 175 10,7 501,54 9,1 2,9
4 Autres céréales 88 5,4 428,99 7,8 4,9
6 Tournesol 76 4,7 319,70 5,8 4,2
16 Fourrage 39 2,4 141,23 2,6 3,6
5 Colza 31 1,9 133,49 2,4 4,3
7 Autres oléagineux 13 0,8 74,82 1,4 5,8
11 Gel (surfaces gelées sans production) 86 5,3 73,82 1,3 0,9
25 Légumes ou fleurs 22 1,4 70,15 1,3 3,2
28 Divers 180 11,1 65,18 1,2 0,4
3 Orge 20 1,2 52,31 1,0 2,6
15 Légumineuses à grains 9 0,6 30,65 0,6 3,4
8 Protéagineux 3 0,2 19,63 0,4 6,5
Total - 1 628 100,1 5 486,78 100,0 3,4
Code
# rm(t1)

1.6 Comparaison avec la répartition des cultures au niveau départemental et national

Code
# jointure spatiale avec la couche département pour récupérer le département où tombe le point

# ouvrir la couche des communes (à convertir en Lambert II epsg 2154)
dep <- s3read_using(
    FUN = sf::read_sf,
    layer = "departement",
    object = "2023/sujet2/diffusion/ign/adminexpress_cog_simpl_000_2023.gpkg",
    bucket = "projet-funathon",
    opts = list("region" = "")) %>% 
  st_transform(2154)

# jointure
df<-point %>% st_join(dep) %>% st_drop_geometry() %>% select(insee_dep)
dep_pt<-df[1,1]
rm(df)

# sinon sélection du département regroupant la plus grande surface agricole
# cas où le cercle recouvre plusieurs départements
# df<-parc_prox %>% st_drop_geometry() %>% 
#   count(insee_dep,wt=surf_parc) %>% 
#   arrange(desc(n))
# dep_pt<-df[1,1]
# rm(df)

# calcul des % surfaces autour du point
stat_pt <- parc_prox %>% st_drop_geometry() %>% 
  count(code_group,wt=surf_parc) %>% add_tally(n) %>% 
  mutate(pct_surf_local=round(100*n/nn,1)) %>%
  select(code_group, pct_surf_local) 
Storing counts in `nn`, as `n` already present in input
ℹ Use `name = "new_name"` to pick a new name.
Code
# récup des % surfaces départementales
stat_dep_pt<-s3read_using(
  FUN=readr::read_rds,
  object = "2023/sujet2/diffusion/resultats/stat_group_cult_by_dep.rds",
  bucket = "projet-funathon",
  opts = list("region" = "")) %>% 
  filter(insee_dep %in% dep_pt) %>% 
  select(insee_dep,code_group,libelle_groupe_culture,pct_surf) %>% 
  rename(pct_surf_dep = pct_surf)

# récup des % surfaces nationales
stat_fm<-s3read_using(
  FUN=readr::read_csv,
  object = "2023/sujet2/diffusion/resultats/stat_group_cult_fm.csv",
  col_types = cols(code_group = col_character()),
  bucket = "projet-funathon",
  opts = list("region" = "")) %>% 
 select(code_group,libelle_groupe_culture,pct_surf) %>% 
  rename(pct_surf_fm = pct_surf)

# appariement des stas locales, départementales, nationales
stat_compar<-stat_fm %>% 
  left_join(stat_dep_pt %>% select(code_group, pct_surf_dep),by="code_group") %>% 
  left_join(stat_pt,by="code_group") %>% 
  select(libelle_groupe_culture,pct_surf_local,pct_surf_dep,pct_surf_fm) %>% 
  arrange(desc(pct_surf_local)) %>% adorn_totals() 

stat_compar %>% 
  setNames(c("Groupe de cultures","surf. locales (%)","surf. départ. (%)","surf. France m. (%)")) %>% kable(
    format="html",
    caption="<span style='font-size:medium'>Comparaison des surfaces locales, départemenales et nationales</span>",
    format.args = list(decimal.mark = ",", big.mark = " "),
    booktabs = TRUE) %>%
  kable_styling(font_size = 15) %>% 
  gsub("font-size: initial !important;",
       "font-size: 20pt !important;",.)%>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>% 
  row_spec(nrow(stat_compar), bold = T, color = "white", background = "grey")
Structure des cultures
Groupe de cultures surf. locales (%) surf. départ. (%) surf. France m. (%)
Maïs grain et ensilage 30,4 19,4 10,0
Prairies permanentes 21,3 29,0 27,7
Blé tendre 13,4 16,4 17,7
Prairies temporaires 9,1 9,5 5,3
Autres céréales 7,8 8,1 3,9
Tournesol 5,8 3,7 2,5
Fourrage 2,6 4,1 3,7
Colza 2,4 2,3 3,5
Autres oléagineux 1,4 0,4 0,7
Gel (surfaces gelées sans production) 1,3 0,6 1,6
Légumes ou fleurs 1,3 1,2 1,5
Divers 1,2 1,0 1,3
Orge 1,0 2,3 6,2
Légumineuses à grains 0,6 0,4 0,2
Protéagineux 0,4 0,9 1,2
Plantes à fibres NA 0,3 0,4
Riz NA NA 0,0
Estives et landes NA 0,1 8,1
Vergers NA 0,1 0,4
Vignes NA 0,2 2,2
Fruits à coque NA 0,0 0,2
Oliviers NA NA 0,0
Autres cultures industrielles NA 0,1 1,8
Total 100,0 100,1 100,1

1.7 Graphique de comparaison des cultures au niveau local, départemental et national

Code
# je sélectionne les 10 groupes de cultures les plus répandus au niveau local 
tab<-stat_compar %>% filter(libelle_groupe_culture!="Total") %>% slice_head(n=10) %>% 
  rename(local=pct_surf_local, departement=pct_surf_dep, france=pct_surf_fm)

# je transpose la table pour rassembler toutes les valeurs dans une seule variable value (ggplot oblige)
tab_piv<-tab %>% pivot_longer(!libelle_groupe_culture) %>% rename(secteur=name) 

# je mets à 0 les valeurs manquantes
#tab_piv<-tab_piv %>% coalesce(value,0L)  
tab_piv[is.na(tab_piv)] <- 0

# levels(as.factor(tab_piv$secteur))
# tab_piv$secteur<-factor(tab_piv$secteur, levels=c("france","departement","local"))
# tab_piv<-tab_piv %>% arrange(desc(secteur), desc(value))

# ggplot => problème de tri : affichage des prairies avant le maïs ?    
p<-ggplot(tab_piv, aes(x = reorder(libelle_groupe_culture,+value), 
                       y = value, 
                       fill=factor(secteur, levels=c("france","departement","local")))) + 
  geom_col(position = "dodge") +
  labs(title="Surfaces comparées des 10 principales cultures locales, en %", x="Culture", y = "%", fill = "Secteur") +
  theme_classic()

# je flippe le graphique pour avoir des barres horizontales  
p+coord_flip()

différences de structure locale, departementale, nationale

1.8 Je teste un graphique par secteur avec facet_wrap

Code
# je sélectionne les 10 groupes de cultures les plus répandus au niveau local 
tab<-stat_compar %>% filter(libelle_groupe_culture!="Total") %>% slice_head(n=10) %>% 
  rename(local=pct_surf_local, departement=pct_surf_dep, france=pct_surf_fm)

# je transpose la table pour rassembler toutes les valeurs dans une seule variable value (ggplot oblige)
tab_piv<-tab %>% pivot_longer(!libelle_groupe_culture) %>% rename(secteur=name) 

# je mets à 0 les valeurs manquantes
#tab_piv<-tab_piv %>% coalesce(value,0L)  
tab_piv[is.na(tab_piv)] <- 0

# ggplot => problème de tri : affichage des prairies avant le maïs ?    
ggplot(tab_piv, aes(x = reorder(libelle_groupe_culture,+value), 
                       y = value)) + 
  geom_col(fill = "lightblue", colour = "black",position = "dodge") +
  labs(title="Surface par culture", x="Culture", y = "%", fill = "Secteur") +
  geom_text(aes(label = value), hjust = -0.3, size = 8/.pt, colour = "black") +
  theme_classic() + coord_flip() + 
  facet_wrap(~secteur,nrow=3,scales='free')

différences de structure locale, departementale, nationale

1.9 Pour aller plus loin, comparer la taille moyenne des parcelles selon les secteurs